% is checked
function y = func_total_nonmod(q, zeta, m, CONSTS)

    check_asympt = false;

    k0 = CONSTS.k0;
    a = CONSTS.a;
    eps_a = CONSTS.eps_a;
    eps = CONSTS.eps;
    eta = CONSTS.eta;
    c = CONSTS.c;
    g = CONSTS.g;
    p  = func_p(q, CONSTS);                
    [q1 q2] = func_qk_from_p_inv(p, CONSTS);
    [delta_1 delta_2] = func_delta_m_l(q, m, CONSTS);
    [alfa_1 alfa_2 beta_1 beta_2] = func_alfa_beta_from_p(p, CONSTS);
    [B_m1 B_m2 B_m1_1 B_m2_1] = func_B_mk_l(q, m, CONSTS);
    [J_m1 J_m2 J_m1_1 J_m2_1] = func_J_mk(q, m, CONSTS);
    [H_m1 H_m2] = func_H_m_l(q, m, CONSTS);
    Q1 = k0*a*q1;
    Q2 = k0*a*q2;
    
    bessel_m_q1 = besselj(m, Q1);
    bessel_mp1_q1 = besselj(m+1, Q1);
    bessel_m_q2 = besselj(m, Q2);
    bessel_mp1_q2 = besselj(m+1, Q2);
              
%     y = (1i*k0/(2*pi)) * (q./(sqrt_from_negative_num(eps_a-q.^2))) .* (exp(-1i*k0*p.*abs(zeta))) .* ...
%         (                                        ...
%           ((1./delta_2) .* ((B_m1 .*                 ...
%              (bessel_mp1_q1 + ...
%                  ((alfa_1.*m.*bessel_m_q1)./(k0*a*q1)))) + ...
%                           (B_m2 .*                 ...
%               (bessel_mp1_q2 + ...
%                  ((alfa_2.*m.*bessel_m_q2)./(k0*a*q2)))))) + ...
%           ((1./delta_1) .* ((B_m1_1 .*                 ...
%              (bessel_mp1_q1 + ...
%                  ((alfa_1.*m.*bessel_m_q1)./(k0*a*q1)))) + ...
%                           (B_m2_1 .*                 ...
%               (bessel_mp1_q2 + ...
%                  ((alfa_2.*m.*bessel_m_q2)./(k0*a*q2)))))) ...
%            );

    y = (1i*k0/(2*pi)) * (q./(sqrt_from_negative_num(eps_a-q.^2))) .* (exp(-1i*k0*p.*abs(zeta))) .* ...
        (                                        ...
          ((1./delta_2) .* ((B_m1 .* J_m1) + (B_m2 .* J_m2_1))) + ...
          ((1./delta_1) .* ((B_m1_1 .* J_m1) + (B_m2_1 .* J_m2_1))) ...
           );
     
%     y1 = (1i*k0/(2*pi)) * (q./(sqrt_from_negative_num(eps_a-q.^2))) .* (exp(-1i*k0*p.*abs(zeta))) .* ...
%           (                                                 ...
%       ((1./delta_2) .* (B_m1 .*                             ...
%              (bessel_mp1_q1 + ...
%                  ((alfa_1.*m.*bessel_m_q1)./(k0*a*q1))))) + ...
%           ((1./delta_1) .* (B_m1_1 .*                       ...
%              (bessel_mp1_q1 + ...
%                  ((alfa_1.*m.*bessel_m_q1)./(k0*a*q1)))))   ...
%                  ); 
             
    y1 = (1i*k0/(2*pi)) * (q./(sqrt_from_negative_num(eps_a-q.^2))) .* (exp(-1i*k0*p.*abs(zeta))) .* ...
          (                                                 ...
      ((1./delta_2) .* (B_m1 .* J_m1)) + ...
          ((1./delta_1) .* (B_m1_1 .* J_m1))   ...
                 ); 
             
%     y2 = (1i*k0/(2*pi)) * (q./(sqrt_from_negative_num(eps_a-q.^2))) .* (exp(-1i*k0*p.*abs(zeta))) .* ...
%           (                                                 ...
%       ((1./delta_2) .* (B_m2 .*                             ...
%              (bessel_mp1_q2 + ...
%                  ((alfa_2.*m.*bessel_m_q2)./(k0*a*q2))))) + ...
%           ((1./delta_1) .* (B_m2_1 .*                       ...
%              (bessel_mp1_q2 + ...
%                  ((alfa_2.*m.*bessel_m_q2)./(k0*a*q2)))))   ...
%                  ); 

    y2 = (1i*k0/(2*pi)) * (q./(sqrt_from_negative_num(eps_a-q.^2))) .* (exp(-1i*k0*p.*abs(zeta))) .* ...
          (                                                 ...
      ((1./delta_2) .* (B_m2 .* J_m2_1)) + ...
          ((1./delta_1) .* (B_m2_1 .* J_m2_1))   ...
                 ); 

    if(check_asympt) % verification of asymptotics - true after q = 150
        
        %---------- alfa_1_s-----------------------------------------------
        figure; plot(q, real(alfa_1), 'b-', q, imag(alfa_1), 'r-'); % asymptotics is hold -> -1
        legend('Re(\alpha_1)', 'Im(\alpha_1)');
        
        %-----------alfa_2_s-----------------------------------------------
        alfa_2_s = -(1/g)*(1-(eta/eps))*q.^2; 
        figure; plot(q, real(alfa_2), 'b.-', q, real(alfa_2_s), 'r-'); title('Re(\alpha_2)'); % asymptotics is hold
        figure; plot(q, imag(alfa_2), 'b.-', q, imag(alfa_2_s), 'r-'); title('Im(\alpha_2)'); % asymptotics is hold (small values)
        
        %-----------beta_1_s-----------------------------------------------
        beta_1_s = ((eta-eps)/(g*eta))*q.^2;
        figure; plot(q, real(beta_1), 'b.-', q, real(beta_1_s), 'r-'); title('Re(\beta_1)'); % asymptotics is hold
        figure; plot(q, imag(beta_1), 'b.-', q, imag(beta_1_s), 'r-'); title('Im(\beta_1)'); % asymptotics is hold (small values)

        %-----------beta_2_s-----------------------------------------------
        beta_2_s = 1 + (g/(eta-eps));
        figure; plot(q, real(beta_2), 'b.-', q, real(beta_2_s), 'r-'); title('Re(\beta_2)'); % asymptotics is hold
        figure; plot(q, imag(beta_2), 'b.-', q, imag(beta_2_s), 'r-'); title('Im(\beta_2)'); % asymptotics is hold (small values)
        
        %----------q1------------------------------------------------------
        q1_s = q;
        figure; plot(q, real(q1), 'b.-', q, real(q1_s), 'r-'); title('Re(q_1)'); % asymptotics is hold
        figure; plot(q, imag(q1), 'b.-', q, imag(q1_s), 'r-'); title('Im(q_1)'); % asymptotics is hold
        
        %----------q2------------------------------------------------------
        q2_s = 1i*(sqrt(-eta/eps))*q;
        figure; plot(q, real(q2), 'b.-', q, real(q2_s), 'r-'); title('Re(q_2)'); % asymptotics is hold (small values)
        figure; plot(q, imag(q2), 'b.-', q, imag(q2_s), 'r-'); title('Im(q_2)'); % asymptotics is hold
        figure; plot(q, abs(q1), 'b.-', q, abs(q2), 'r-'); title('|q_1| & |q_2|');

        %----------p-------------------------------------------------------
        p_s = -1i*q;
        figure; plot(q, real(p), 'b.-', q, real(p_s), 'r-'); title('Re(p)'); % asymptotics is hold
        figure; plot(q, imag(p), 'b.-', q, imag(p_s), 'r-'); title('Im(p)'); % asymptotics is hold
        
        %----------n1------------------------------------------------------
        [n1 n2] = func_nk_from_p(p, CONSTS);
        n1_s = -(g./p).*(eta/(eta-eps));
        figure; plot(q, real(n1), 'b.-', q, real(n1_s), 'r-'); title('Re(n_1)'); % asymptotics is hold
        figure; plot(q, imag(n1), 'b.-', q, imag(n1_s), 'r-'); title('Im(n_1)'); % asymptotics is hold
        
        %----------n2------------------------------------------------------
        n2_s = ((eta-eps)/g)*p;
        figure; plot(q, real(n2), 'b.-', q, real(n2_s), 'r-'); title('Re(n_2)'); % asymptotics is hold (small values)
        figure; plot(q, imag(n2), 'b.-', q, imag(n2_s), 'r-'); title('Im(n_2)'); % asymptotics is hold
        
        %----------(alpha_1+beta_1)/Q1^2-----------------------------------
        expr1 = (alfa_1 + beta_1)./((k0*a*q1).^2);
        expr1_s = (1/((k0*a)^2))*((eta-eps)/(eta*g));
        figure; plot(q, real(expr1), 'b.-', q, real(expr1_s), 'r-'); title('Re((\alpha_1 + \beta_1)/Q_1^2)'); % asymptotics is hold
        figure; plot(q, imag(expr1), 'b.-', q, imag(expr1_s), 'r-'); title('Im((\alpha_1 + \beta_1)/Q_1^2)'); % asymptotics is hold

        %----------n_1*n_2-------------------------------------------------
        expr2 = n1.*n2;
        expr2_s = -eta;
        figure; plot(q, real(expr2), 'b.-', q, real(expr2_s), 'r-'); title('Re(n_1*n_2)'); % asymptotics is hold
        figure; plot(q, imag(expr2), 'b.-', q, imag(expr2_s), 'r-'); title('Im(n_1*n_2)'); % asymptotics is hold (small values)

        %----------Delta_1-------------------------------------------------
        delta_1_s = 1i*((eta-eps)/g)*(eta/eps_a)*(1./(((k0*a)^2)*(sqrt(eta/eps))*q)).*...
        ((bessel_mp1_q1./bessel_m_q1)+1i).*((bessel_mp1_q2./bessel_m_q2)+...
        (eps_a/eta)*(sqrt(eta/eps))*1i);
        
        figure; plot(q, real(delta_1), 'b.-', q, real(delta_1_s), 'r-'); title('Re(\Delta^{(1)})'); % numerical function has the singular peaks
        figure; plot(q, imag(delta_1), 'b.-', q, imag(delta_1_s), 'r-'); title('Im(\Delta^{(1)})'); % asymptotics is hold

        %----------Delta_2-------------------------------------------------
        delta_2_s = -1i*((eta-eps)/g)*(eta/eps_a)*(1./(((k0*a)^2)*(sqrt(eta/eps))*q)).*...
        ((bessel_mp1_q1./bessel_m_q1)-1i).*((bessel_mp1_q2./bessel_m_q2)-...
        (eps_a/eta)*(sqrt(eta/eps))*1i);
        
        figure; plot(q, real(delta_2), 'b.-', q, real(delta_2_s), 'r-'); title('Re(\Delta^{(2)})'); % numerical function has the singular peaks
        figure; plot(q, imag(delta_2), 'b.-', q, imag(delta_2_s), 'r-'); title('Im(\Delta^{(2)})'); % asymptotics is hold
        
        %----------Bm_1----------------------------------------------------
        B_m1_s = ((4*pi)/c)*(k0*a).*...
                 (n2.*(eta/eps_a).*(J_m2-(eps_a/eta)*H_m2).*H_m2);
        figure; plot(q, real(B_m1), 'b.-', q, real(B_m1_s), 'r-'); title('Re(B_{1m})'); % asymptotics is hold
        figure; plot(q, imag(B_m1), 'b.-', q, imag(B_m1_s), 'r-'); title('Im(B_{1m})'); % asymptotics is hold
        
        %----------Bm_2----------------------------------------------------
        B_m2_s = ((4*pi)/c)*((k0*a)).*...
                  1i.*(eta/eps_a).*(m./(((k0*a)^2)*q)).*(J_m1-H_m2);
        figure; plot(q, real(B_m2), 'b.-', q, real(B_m2_s), 'r-'); title('Re(B_{2m})'); % numerical function has the singular peaks
        figure; plot(q, imag(B_m2), 'b.-', q, imag(B_m2_s), 'r-'); title('Im(B_{2m})'); % asymptotics is hold

        %----------Bm_1_1--------------------------------------------------
        B_m1_1_s = ((4*pi)/c)*((k0*a)).*...
                    (n2.*(eta/eps_a).*(J_m2-(eps_a/eta)*H_m1).*H_m1);
        figure; plot(q, real(B_m1_1), 'b.-', q, real(B_m1_1_s), 'r-'); title('Re(B_{1m}^{(1)})'); % asymptotics is hold
        figure; plot(q, imag(B_m1_1), 'b.-', q, imag(B_m1_1_s), 'r-'); title('Im(B_{1m}^{(1)})'); % asymptotics is hold
        
        %----------Bm_2_1--------------------------------------------------
        B_m2_1_s = ((4*pi)/c)*((k0*a)).*...
              1i.*(eta/eps_a).*(m./(((k0*a)^2)*q)).*(J_m1-H_m1);
        figure; plot(q, real(B_m2_1), 'b.-', q, real(B_m2_1_s), 'r-'); title('Re(B_{2m}^{(1)})'); % numerical function has the singular peaks
        figure; plot(q, imag(B_m2_1), 'b.-', q, imag(B_m2_1_s), 'r-'); title('Im(B_{2m}^{(1)})'); % asymptotics is hold

        %----------y1------------------------------------------------------
        y1_s = func_1_nonmod_analytical(q, zeta, m, CONSTS);
        figure; plot(q, real(y1), 'b.-', q, real(y1_s), 'r-'); title('Re(the first summand)'); % asymptotics is hold (small values)
        figure; plot(q, imag(y1), 'b.-', q, imag(y1_s), 'r-'); title('Im(the first summand)'); % asymptotics is hold
    
        %----------y2------------------------------------------------------
        y2_s = func_2_nonmod_analytical(q, zeta, m, CONSTS);
        figure; plot(q, real(y2), 'b.-', q, real(y2_s), 'r-'); title('Re(the second summand)'); % asymptotics is hold (small values)
        figure; plot(q, imag(y2), 'b.-', q, imag(y2_s), 'r-'); title('Im(the second summand)'); % a little occilations for full func
                                                                                                % an exp decrising for approximation
    end
       
end
      
        







